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Abstract 


This  paper  describes  a  method  of  simulating  a  Markov  chain  for 
the  purpose  of  estimating  functions  of  the  chain  and  functions  of 
associated  semi-Markov  processes.  In  particular,  special  attention  is 
devoted  to  the  estimation  of  the  probability  density  function  of  first 
passage  time  from,  say,  state  a  to  state  b.  Rotation  sampling  Is  used 
to  achieve  variances  of  estimators  of  order  where  k  is  the 

number  of  replications,  which  compares  with  0(l/k)  when  Independently 
sampled  replications  are  used.  Since  both  Independent  and  rotation  sampling 
have  computation  time  complexity  0(k),  the  relative  advantage  of  rotation 
sampling  Is  clear  as  k  -*■  «.  The  paper  presents  two  examples  to  Illustrate 


the  method. 


•td 


Consider  a  positive  recurrent  aperiodic  Markov  chain  with  state  space 
S  =  (0,1,..., n)  and  transition  probability  matrix  =  ||  p^^  ||  where 
p^j  denotes  the  probability  of  moving  from  state  i  to  state  j  in  one 
step  for  i,j  €  S  .  Let  s..  denote  the  frequency  of  one-step  transitions 

'  J 

from  i  to  j  during  a  first  passage  from  state  a  to  state  b  a,b  e  S 
and  let  s  denote  a  lx(n+l)  vector  with  s . .  in  column  n(i-l)  +  j  . 

•  J 

Also  let  p(s;a,b)  denote  the  probability  of  observing  the  frequency 
vector  s  during  a  first  passage. 

The  need  to  know  p(s;a,b)  arises  in  many  areas  of  applied  probabilitv 
and  operations  research.  For  example,  consider  a  semi-Markov  process  with 
Markov  chain  and  holding  time  probability  density  function  (p.d.f.) 
f^j(t)  on  [0,-)  for  one-step  transitions  from  i  to  j  i , j  c  S  .  Let 
f..(t|s)  denote  the  s-fold  convolution  of  fi<(t)  with  itself.  Then  the 

I  J  I J 

semi-Markov  process  has  first  passage  time  p.d.f. 

g(t;a,b)  =  p(£;a,b)  f(t|j^)  (1) 

where  f(t|£)  denotes  the  convolution  of  the  p.d.fs.  f.^(t|s..)  for  all 

I  J  I  J 

i,j  €  S  such  that  s. ,  >  0  . 

Although  one  can  write  down  systems  of  equations  to  represent  (1),  it  is 
difficult  to  solve  these  systems.  In  particular,  it  is  difficult  to  derive 
computationally  convenient  expressions  for  p(s;a,b)  .  The  purpose  of 
this  paper  is  to  propose  a  method  for  estimating  p(s;a,b)  by  means 
of  a  Bupereffiaient  Monte  Carlo  sampling  method  and  then  showing  how  these 
estimates  can  be  used  to  estimate  an  important  class  of  reward  functions,  of 
which  (1)  is  a  special  case.  In  k  independent  Monte  Carlo  replications  or 
sampling  experiments,  the  variances  of  estimators  usually  are  0(1 /k).  The 
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sampling  plan  proposed  here  leads  to  variances  of  order  0(l/k  ).  Since  the 
computation  time  complexity  of  the  proposed  scheme  is  0(k),  which  also  holds 
for  independent  replications,  the  proposed  sampling  plan  is  termed  supeveffiaient 
The  sampling  plans  use  a  special  case  of  the  antithetic  variate  method 
(Hammersley  and  Handscomb  1964)  called  rotation  sampling  which  is  described 
in  detail  in  Fishman  and  Huang  (1980).  The  technique  has  already  been 
applied  to  the  direct  simulation  of  Markov  chains  in  Fishman  (1981a,  1981b, 
1982a).  The  present  paper  extends  the  technique  to  the  estimation  of  the 
more  complex  probabilities  p(s;a,b). 

Section  1  describes  the  procedure  in  detail  and  Section  2  describes 
how  it  applies  for  Gamma  distributed  holding  times.  Section  3  presents  two 
examples,  one  of  which  is  the  estimation  of  the  p.d.f.  of  busy  time  for  the 
M/M/1  queueing  model.  The  results  of  this  example  are  compared  with  the 
theoretical  solution  in  Cox  and  Smith  (1961,  p.  148).  These  examples  use 
the  MCHAIN  FORTRAN  subprogram  (Fishman  1982b)  which  enables  one  to  simulate 
k  replications  of  the  chain  in  parallel,  using  rotation  sampling. 

For  consistency  with  the  earlier  work  in  Fishman  (1981a,  1981b,  1982a), 
we  add  several  modelling  details.  We  assume  that  there  exists  a  positive 
integer  «  <  (n-l)/2  such  that 


and 


Pij 


for  li-j  I  >  6 


rmin(n,i+6)  _  , 

^j=max(0,i-6)  ^ij  ”  ‘ 

Also,  let  Sj  denote  the  total  number  of  states  that  have  positive  transition 
probabilities  from  state  j  and  let  r  =  l.*-**Sj}  denote  the  ordered 

sequence  (m.  <  m.  r  =  1 . s.  -1)  of  the  s.  states  to  which  entry  can 

J  *  J  •  *  ^  »  J  J 

occur  from  state  j.  Then  one  has  the  representation 
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r  “  1 1 . . .  ,s j 


I  ^  P-  =  1 
J"'jr 


6  ^  max  (| mj  j|  .  |  m^.^  -j  | ) 

3 


j=  0,1 . 


The  value  of  this  alternative,  but  equivalent,  representation 
becomes  apparent  when  actually  generating  sample  paths  by  simulation  on 
a  computer.  See  Fishman  (1982b). 

1 .  Simulation  of  a  Markov  Chain 

Reward  Functions 

Let  I  denote  the  set  of  all  one-step  (i.j)  transition  pairs  with  p..>  0. 

^  J 

Suppose  one  wants  to  use  simulation  to  study  the  behavior  of  the  chain 

during  an  interval  that  begins  with  exit  from  state  a  and  ends  with  entry  to 

state  b  where  a,b  e  S.  Let 

N.  =  frequency  of  transitions  of  type  i  during 
^  first  passage  on  an  arbitrarily  selected  (3a) 
sample  run 
and 

A(j,,...,j  ;a)  =  reward  received  during  first  passage  for 

transitions  of  types  l,...,r  (3b) 

respectively. 

A  type  denotes  all  the  one-step  (i,J)  transition  pairs  which  lead 
to  the  exact  same  reward.  Here  there  are  r  distinct  types.  This  notation 
proves  convenient  later  in  reducing  the  dimensionality  of  the  estimation 
problem. 


On  k  replications  the  sample  mean  reward  is 

R  =  —  A(N^'^^ 


where  the  superscript  (£)  denotes  replication  t.  The  reward  function  A  can 
assume  many  alternative  forms.  For  example,  if 

A(j^,...,j^;£)  »  j^a^  (5) 
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for  given  S,  *  (a^,...,a^)  ,  Rj^  is  linear  in  where 

Sj  =  n/>.  (6) 

This  case  has  been  studied  in  detail  from  the  point  of  view  of  simulation  in 
Fishman  (1981a,  1981b,  1982a).  As  a  second  example,  consider  a  semi -Markov 
process  with  Markov  chain  and  continuous  holding  time  p.d.fs. 

^]>*--»f^  the  r  types  of  transitions  in  I.  Let 


g^(tlj)  =  f^.(t)  *  f^(t)  *  ...  *  f^(t)  0  £  t  <  »  (7) 

denote  the  j-fold  convolution  (j  >  0)  of  f^.  .  Then  for  given 
i  =  l,...,r  the  p.d.f.  of  first  passage  time  from  a  to  b  is 

A(j^,...,j^;0)  =  g^(tlj^)  *  g2(t|j2)  *...*  g^(tlj^)  (8) 

where  the  convolution  includes  only  those  g.'s  for  which  j.  >  0  and  e  =  t 

Aggregation  Across  Replications 

Although  (4)  proves  useful  for  computation  of  Rj^,  an  alternative 
representation  considerably  simplifies  the  derivation  of  results.  Let 


■J, . "s=l' 


where 


6(x)  =  1 
=  0 


if  X  =  0 
otherwise. 


so  that  one  can  write  (4)  as 


Note  that 


1  00 

'  k  . 

>  _  rm  1. 
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is  the  number  of  replications  that  enter  b  on  step  m.  More  importantly,  note 
that  K.  .  /k  is  an  unbiased  estimator  of  p(s;a,b)  in  (1)  with 

J] . Jy.  ~ 

frequency  vector  s  =  in  the  case  in  which  there  are  r  =  (n  +  1)^ 

distinct  types  of  transitions.  Although  we  hereafter  discuss  the  properties 
of  R|^  it  should  be  recognized  that  a  reward  function 

A(ji , . . .  "  j\j) 

specializes  the  analysis  to  the  estimation  of  p(s;a,b)  where  s  = 

Serial  Simulation 

Presumably,  the  objective  of  simulation  is  to  perform  k  replications 
sufficient  to  achieve  an  acceptable  accuracy  for  R|^  as  an  estimate  of  ERj^. 

In  the  case  of  k  independent  replications  in  series,  with  a^b,  one  can  ensure 
that  each  replication  begins  with  an  exit  from  a  and  ends  with  a  first  entry 

into  b  by  replacing  J*0,1 . n)  by  and  for  V  Sj^a. 

For  the  case  a*b  there  is  no  need  to  modify  the  probabilities.  In  the 
independent  case  one  has  var  R|^=0(l/k)  and  computation  time  complexity  0(k). 
Parallel  Simulation 

To  speed  up  the  convergence  rate  of  var  Rj^  we  turn  to  parallel 

simulation  using  rotation  sampling.  Hereafter,  a  prime  superscript  denotes 

parallel  rotation  sampling.  Consider  k  replications  each  beginning  with  an 

exit  from  state  a  and  set  Pjjj,*!  and  Pjjj'O  for  i  sj<b.  If  a»b  this  modification 

of  tPjjjl  should  be  performed  after  the  first  transition.  This  modification 

makes  the  chain  absorbing  with  transient  states  where  =  S-b, 

and  absorbing  state  b.  Let  T„  denote  the  set  of  transient  states  with  at 

m 

least  one  resident  replication  at  the  end  of  step  m  and  let  K.'  denote  the 

jm 

number  of  such  replications  in  state  jcT  .  Let  ;  JeT  )  denote  a  sequence 

m  jm  m  ^ 

of  i.i.d.  random  variables  uniformly  distributed  on  [0,1)  and  define  for 
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U.  <l-(i-l)/K. 
jm  jtn 


=  U.  +  (i-l)/K.' -1  otherwise, 
jm  '  jm 


(10) 


Also  define 


and 


‘’js  "  ^1=0  Pj£ 


scS 


(11) 


Then  at  step  m+1  replication  i  in  state  j  goes  to  state  s  if 

q.  ,  <  V..  <  q.  .  Since  q,.  =1,  a  destination  is  guaranteed. 

J.s-1  -  TOtn  js  ^jn 

Note  that  althouqh  V,.  ....,Vw.  are  each  unlfomly  distributed  on 

ijm 

[0,1),  they  are  not  independent.  Therefore,  the  sample  path  of  each  replication 
from  a  to  b  has  the  correct  probability  law,  but  the  k  paths  are  not  independent. 
The  assignment  in  (10)  is  called  rotation  sampling  and  produces  a  clear  benefit 
in  simulation.  For  a  reward  function  as  in  (5),  the  use  of  (4)  leads  to 
var  _<  0((ln  k/k)^)  n  <  » 

and  (12) 

var  _<  0((ln  k)^/k^)  n  »  , 

which  Improve  on  0(1 /k)  for  Independent  replications.  Moreover, 
for  (5)  one  can  show  that  the  computation  time  complexity  is  0(ln  k)  for  n  <  » 
and  0((ln  k)  )  for  n  -»■  «*.  This  last  result  follows  from  a  more  compact 
sampling  scheme,  than  that  in  (10),  that  can  be  used  when  there  is  no  need 
to  keep  track  of  the  distinct  sample  paths  on  each  replication.  See 
Fishman  (1981a,  1981b,  1982a).  For  the  more  general  reward  function  (4), 
keeping  sample  path  data  to  compute  jel)  for  k  is  necessary  and 

therefore  the  computation  time  complexity  is  0(k)  . 


One  additional  result  is  of  importance.  Let  Kl.  denote  the  number  of 

1  jm 

replications  that  move  from  i  to  j  on  step  m.  Then  it  is  shown  in  Fishman  (1981a) 
that  as  !(-»•“> 

var  =0(1)  iJcS  m=l,2 .  (13) 

We  use  this  result  next  when  bounding  var  Rj^  for  the  general  reward  function  (4). 
Observe  that  by  analogy  with  (9) 


so  that  for  given  m  and  k  ^  « 


Zv....i,=o  ZV....VO  . 


• » j  ) 


il+...+i^=m  ji+...+j^*m 


=  0(1)  . 


Then  i  t  can  be  shown  that  as  k  • 


cov  (Kj  i  •  Ki  .  )  =  0(1)  (16) 

so  that  our  estimate  Kl  K\  /k  of  p(s;a,b)  with  s  *  (j, . j  )  hi 

J] » •  ♦  •  >  ~  I  r 

variance  0(l/k^)  .  Proof  of  (16)  follows  from  Proposition  1. 


Proposition  1. 


Let  X^,...,X  denote  random  variables  with 
var  X^  »  0^(2)  , 
cov  (X^.  Xj)  =  o^j(z) 


and  define 


with 


If 


then 


Proof.  Observe  that 


where 


2 

var  Y.  =  uiAz)  . 

J  J 

2 

litn  0).  (z)  =  c  <  »  V  j, 
z-*®  J 


2 

litn  cr.(z)  ^constant  Vi.  (17) 

Z-vo  ^ 


a^(z)  =  “t^z)  -  u^_i(z)  -  2C^(z) 


C^(z)  =  cov  (X^.  t>l  . 


Since 

one  has 

o^{z)  -  2  o^(z)  (*)^_^(z)  +  u^_-j  (z)  -  a)^(k)  * 

[at(z)  -  w^_i{z)  -  a)^(z)]  [o^(z)  -  (*>^_i(z)+  u)^(z)]  <  0 

SO  that 

«t,i(2)  -  -  °t^^^  -  “t-1^^^  ■*■  "t^^^ 

Taking  the  limit  as  z-*- <»  gives  (17).  The  result  in  (16)  follows  directly. 
We  now  derive  var  Rj^  for  the  important  case 

. j^=0  S)  <  constant  (19a) 

and 

A(ji»... » Jy,  !^)  V  j'|f...>J^. 


(19b) 
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Suppose  that  a  transition  of  type  iei  has  the  Gamma  p.d.f. 


r  (a^.) 


a.  >  1  ,  X.  >  0, 
1  ~  1 

0  <  t  <  * 


(24) 


denoted  by  G(a^,  x^. )•  Let  us  also  assume  that  all  transition  pairs  {i,j) 
with  the  same  p.d.f.  are  regarded  as  one  type  in  I  and  that  x^<  ...<  x^. 

For  an  arbitrary  replication  with  transitions  of  type  i,  the  corresponding 

total  holding  time  has  the  Gamma  distribution  6{a.  N^. ,  x^.)  with  p.d.f. 


gi(t|N^)= 


x“i^• 


“i*  >  1 ,  >  0, 

0  <  t  <  »  .  (25) 


Then  the  estimated  first  passage  time  p.d.f.  is 


where 


9(tlj^,...,jr) 


(26a) 


J^)  =  gi(tlj^)  *...*  g^(tjj^) 


(26b) 


and  the  convolution  is  taken  only  over  those  types  i  for  which  >  0  . 

In  practice,  it  is  more  convenient  to  compute  R|[(t)  from 


(27) 


which  is  algebraically  equivalent  to  (26a).  Note  that  the  sample  frequencies  (9) 
for  the  parallel  simulation,  need  not  be  computed  when  (27)  is  used.  However, 
the  computation  of  g(tiN^^'^\...,N^'“^)  remains  a  difficult  one  when 

are  unequal.  By  using  the  Taylor  series  expansion  of  e®  ,  one 
can  write  (27)  equivalently  as 
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ProposHion  2.  Let  the  reward  function  A  satisfy  (19a)  and  (19b).  Then 


var  Rj^  <  0(l/r). 


(20) 


Proof.  Let  B=sup 

>  •  •  •  ij, 

bounded,  clearly 


var  K.  ...  Since  these  variances  are 

J.|  .  »jp 


var  R|^  <  . A(j^ . j^)]^  B/k^ 


(21) 


<  0(l/k^). 

We  now  show  that  the  result  (20)  applies  for  the  case  of  the  first 
passage  time  p.d.f.  in  (8).  Observe  that  one  can  write 

L  ,*  =0  ^(j]  ••  •  •  »"ip* 


h^{t)  *  h2(t)  *...*  h^(t) 


(22) 


where  h^  is  the  renewal  density  function 


h,(t)  •  g,(t(j)  0  i  t  <  -  (23) 

and 

g^.(tlO)=  0  i=l,...,r. 

Since  each  of  the  renewal  densities  can  be  bounded  from  above,  it  is  clear 
that  for  t€  [0,®)  (22)  is  also  bounded,  thus  establishing  (20). 

2.  Gamma  Distributed  Holding  Times 

Although  our  estimator  of  p(5,;a,b)  based  on  rotation  sampling  has 

a  clear  advantage  over  an  estimator  based  on  independent  replications, 
computational  problems  remain  in  using  either  estimator  to  estimate  the 
first  passage  time  p.d.f.  in  (1).  In  this  section  we  describe  and  illustrate 
how  these  problems  can  be  overcome  for  the  common  case  of  Gamma  holding  times. 


where 


^is_-|=0  9s^^'”s  ^s-1^ 


s  =  max  (i:  N.  >  0), 


^0  =  0 


and  for  j  =1 ,. . . ,s 


’j  "  'j-l  *  ’o’ 

"j  •  1^.1 


if  N.=0 

J 


“j  "  A, 


t  -  min(i:  i  >  j,  >  0). 


More  concisely,  (28)  is  equivalent  to 


g(tlN 


where 


=  ^i^_^=0  S-1  ^^r-1^  9s^^l”s  ^s-1^ 


did])  =  Pi(1i  |0) 


Finally  (28)  has  the  equivalent  form 

,  ,  (Ajt)  V*st 

9(t  », . N  •— = - — - q.  ,(0)  t  y*  ^ - 

r  »  ,  •'s-l  ‘’=1^  („  t£.i) 


T 


Expression  (31)  is  exact  and  can  be  used  in  (27)  with  N^.  '  replacing  N^. 

i  =  l,...,r  to  compute  an  estimate  of  the  p.d.f.  Note  that  the  term  in 
brackets  in  (31)  is  eventually  a  decreasing  function  of  j  for  given  t  , 
This  fact  enables  one  to  truncate  the  summation  on  j  at  a  point  where  the 
remainder  is  relatively  incidental. 

3.  Examples 

Consider  an  M/M/1  queueing  model  v/ith  arrival  rate  x  and  service  rate  i 
Suppose  one  wants  to  estimate  the  p.d.f.  of  the  busy  period.  Cox  and 


Smith  (1961,  p.  148)  give  this  p.d.f.  as 


g(t)  * 


t(x/oj) 


I-j  (2t  /xHi) 


0  <  t  <  ® 


where  I-j  denotes  the  Bessel  function  of  imaginary  argument  and  first  order. 

In  this  example,  a=b=0,  n=«  and  r=l,  since  all  holding  time  distributions 

from  states  i>0  are  exponential  with  identical  rate  x+w.  If  N  is  the  total 

number  of  transitions  on  an  arbitrarily  selected  replication,  excluding 

the  initial  transition  from  state  0,  then 

/ , .  \N  aN"!  _“(x+(o)t 

A(N,t)  =  g(tlN)  =  ^ ^ -  (33) 

r(N) 

and  the  sample  mean  based  on  parallel  simulation  is 

-(x+u))t 


K'.  =  6(N^^^-j). 
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Experimental  Design 

We  distinguish  between  mlcAaKtpUjcuxtioYit  and  macAo/iepttca^oni. - 
In  particular,  each  of  L  Independent  macrorepll cations  of  the  simulation 
contains  k  parallel  ml crorepll cations  that  use  (10).  The  L  Independent  data 
blocks  enable  us  to  estimate  var  R^(t)  for  each  value  of  k  considered.  The 
global  factor  levels  are  x=0.9,  o)=l  and  t=0.5j;  j=l,,..,l0.  For  these,  L=100 
macrorepllcatlons  were  run  for  six  experiments  with  k.=2^'*''^  microrepllcatlons 

J 

on  experiments  j=l,...,6  .  The  macrorepllcatlons  were  used  to  estimate 
var  R|^(t)  for  each  t  and  k  considered.  Also,  for  these  x,u  and  t's, 
L=1000  macrorepllcatlons  were  run  each  for  k=l  microrepllcatlon  to  get  a 
baseline  estimate  of  var  Rj(t)  ,  the  variance  without  rotation  sampling. 

Table  1  shows  g(t)  from  (32)  and  ^256^^^  comparative  purposes. 

Insert  Table  1  about  here. 

To  measure  variance  reduction  we  use  the  ratio 

var  R;(t) 

V  s  _ ! _ 

^  k  var  Rj[(t) 

Our  earlier  results  Indicate  that  this  quantity  Is  linear  In  k  as  k-*^  . 

Now  when  assessing  variance  reduction  for  the  estimation  of  a  function,  it 
should  be  noted  that  the  variance  reductions  for  all  10  values  of  t  in 
Table  1  need  not  be  Identical  for  each  k  .  Therefore  for  each  k  we 

plotted  the  minimal  and  maximal  Vj^'s  In  Fig.  1,  noting  that  all  remaining 

variance  reductions  lie  between  them.  In  general,  the  graphs  agree  with 

theory.  Moreover,  we  observed  that  for  fixed  k  the  V. 's  decreased  with 

k 

Insert  Figure  1  about  here. 
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decreasing  coefficients  of  variations  computed  for  the  estimates  derived  in 
the  baseline  case  of  independent  replications.  This  result  is  reassuring 
for  it  indicates  that  the  benefits  of  rotation  sampling  are  most  apparent 
precisely  where  needed,  namely  where  the  larger  coefficients  of  variation 
arise. 

Our  second  example  illustrates  the  estimation  technique  based  on  (31). 

It  uses  the  same  Markov  chain  as  in  example  1  but  now  with  x=0.5,  a)=l 
and  holding  time  distributions  G(l,x^)  m=l,...,4  .  Table  2  lists  the 
transition  types.  The  objective  was  to  estimate  the  first  passage  time 

Insert  Table  2  about  here. 

for  a=b*0  at  times  t»0.1j  j=l,...,40  .  For  these  parameters,  L=100 

macroreplications  were  run  for  nine  experiments  with  k.=2^'’^'^  microreplications 
on  experiment  j=l,...,9  .  Also,  L*1000  macroreplications  were  run  each  with 
k=l  micrereplication  to  get  a  baseline  estimate  of  var  R.j(t)  .  Figure  2 
shows  the  graphs  of  minimal  and  maximal  variance  reduction  V|^  versus  k  on 
logarithmic  scales.  Again  the  graphs  essentially  agree  with  theory.  For  each  k 
the  same  relationship  between  variance  reduction  and  coefficients  of  variation 
were  observed  as  in  example  1. 

4.  Conclusions 

This  paper  demonstrates  how  one  can  estimate  by  a  superefficient 
method  the  probability  of  observing  s^,...,Sj.  transitions  of  types  l,...,r 
during  a  first  passage  from  state  a  to  state  b  in  a  positive  recurrent 
aperiodic  Markov  chain  with  state  space  S  -  (0,1,..., n)  .  It  then  shows  how 
these  results  can  be  used  to  estimate  a  class  of  mean  reward  functions  of  which 


first  passage  time  p.d.fs.  are  a  special  case.  Section  2  specializes  the 
analysis  to  the  case  of  Gainna  holding  times,  showing  how  one  derives  expression 
(31)  which  is  convenient  for  numerical  calculation.  Section  3  demonstrates 
the  technique  by  two  examples. 

In  order  to  implement  the  proposals  of  this  paper,  one  needs  a  package 

It) 

that  generates  the  transition  frequencies  ’  ;  i,j  e  S}  £=!,..., k 

for  k  parallel  replications  based  on  rotation  sampling.  The  MCHAIN  FORTRAN 
program  in  Fishman  (1982b)  will  generate  these  data. 
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